/*
 * Copyright (C) 2008 Peter Kling
 *
 * This file is part of 'A Kind of Billiard', a graphical billiard simulation.
 *
 * 'A Kind of Billiard' is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option) any
 * later version.
 *
 * 'A Kind of Billiard' is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
 * more details.
 *
 * You should have received a copy of the GNU General Public License along with
 * 'A Kind of Billiard'. If not, see <http://www.gnu.org/licenses/>.
 */

package cg2.mathematics;

/**
 * @brief This class collects some global mathematical methods.
 */
public final class Mathematics {
	// -----------------------------------------------------------------------------------------------------------------
	// Distance methods
	// -----------------------------------------------------------------------------------------------------------------
	
	/**
	 * @brief Returns the distance between two points.
	 */
	public static double distance(MathVector p1, MathVector p2) {
		return Mathematics.sub(p2, p1).norm();
	}
	
	/**
	 * @brief Returns the distance between the given point and plane.
	 */
	public static double distance(MathVector p, MathPlane e) {
		return Math.abs(distance_oriented(p, e));
	}
	
	/**
	 * @brief Returns the distance between the given vector and plane.
	 */
	public static double distance(MathPlane e, MathVector v) {
		return distance(v, e);
	}
	
	/**
	 * @brief Returns the oriented distance between the given point and plane.
	 * 
	 * If the oriented distance is > 0, the given point lies in the half space containing the origin.
	 */
	public static double distance_oriented(MathVector p, MathPlane e) {
		return e.getDistance() - Mathematics.mul(p, e.getNormal());
	}
	
	/**
	 * @brief Returns the oriented distance between the given point and plane.
	 * 
	 * If the oriented distance is > 0, the given point lies in the half space containing the origin.
	 */
	public static double distance_oriented(MathPlane e, MathVector p) {
		return distance_oriented(p, e);
	}
	
	// -----------------------------------------------------------------------------------------------------------------
	// Misc operations
	// -----------------------------------------------------------------------------------------------------------------
	
	/**
	 * @brief Returns the reflection of the given point along the given plane.
	 */
	public static MathVector mirrorPoint(MathVector p, MathPlane e) {
		double distance_oriented = e.getDistance() - Mathematics.mul(p, e.getNormal());
		return add(p, mul(e.getNormal(), 2*distance_oriented));
	}
	
	/**
	 * @brief Returns the reflection of the given point along the given plane.
	 */
	public static MathVector mirrorPoint(MathPlane e, MathVector p) {
		return mirrorPoint(p, e);
	}
	
	/**
	 * @brief Returns the reflection of the given vector along the given plane.
	 * 
	 * The reflection of a vector along a plane is defined by the vector resulting from mirroring v as a point along the
	 * translation of e into the origin.
	 */
	public static MathVector mirrorVector(MathVector v, MathPlane e) {
		MathPlane e_prime = new MathPlane(e.getNormal(), 0.);
		return mirrorPoint(v, e_prime);
	}
	
	/**
	 * @brief Returns the reflection of the given vector along the given plane.
	 * 
	 * The reflection of a vector along a plane is defined by the vector resulting from mirroring the point v along the
	 * translation of e into the origin.
	 */
	public static MathVector mirrorVector(MathPlane e, MathVector v) {
		return mirrorVector(v, e);
	}
	
	// -----------------------------------------------------------------------------------------------------------------
	// Vector operations
	// -----------------------------------------------------------------------------------------------------------------

	/**
	 * @brief Returns the sum of two vectors.
	 * @throws IllegalArgumentException if the dimensions of the vectors do not match. 
	 */
	public static MathVector add(MathVector v1, MathVector v2) {
		if (v1.dim != v2.dim)
			throw new IllegalArgumentException("Dimensions of vectors do not match.");
		MathVector v = new MathVector(v1.dim);
		v.add(v1);
		v.add(v2);
		return v;
	}
	
	/**
	 * @brief Returns the i-th n-dimensional unit vector.
	 */
	public static MathVector getUnitVector(int n, int i) {
		MathVector v = new MathVector(n);
		v.set(i-1, 1.);
		return v;
	}
	
	/**
	 * @brief Returns the component wise product of the given vector and scalar.
	 */
	public static MathVector mul(double s, MathVector v) {
		return Mathematics.mul(v, s);
	}

	/**
	 * @brief Returns the component wise product of the given vector and scalar.
	 */
	public static MathVector mul(MathVector v, double s) {
		MathVector result = v.clone();
		result.mul(s);
		return result;
	}

	/**
	 * @brief Return the scalar product of two vectors.
	 * @throws IllegalArgumentException if the dimensions of the vectors do not match.
	 */
	public static double mul(MathVector v1, MathVector v2) {
		if (v1.dim != v2.dim)
			throw new IllegalArgumentException("Dimensions of vectors do not match.");
		double result = 0;
		for (int i=0; i<v1.dim; i++)
			result += v1.get(i)*v2.get(i);
		return result;
	}

	/**
	 * @brief Returns the difference of two vectors.
	 */
	public static MathVector sub(MathVector v1, MathVector v2) {
		if (v1.dim != v2.dim)
			throw new IllegalArgumentException("Dimensions of vectors do not match.");
		MathVector v = new MathVector(v1.dim);
		v.add(v1);
		v.sub(v2);
		return v;
	}
	
	// -----------------------------------------------------------------------------------------------------------------
	// Matrix operations
	// -----------------------------------------------------------------------------------------------------------------
	
	/**
	 * @brief Returns the addition of the given matrices.
	 * @throws IllegalArgumentException if matrix dimensions do not match.
	 */
	public static MathMatrix add(MathMatrix m1, MathMatrix m2) {
		if (m1.n != m2.n || m2.m != m2.m)
			throw new IllegalArgumentException("Given matrices can not be added (dimensions missmatch).");
		MathMatrix result = m1.clone();
		result.add(m2);
		return result;
	}
	
	/**
	 * @brief Returns a rotation matrix for rotations of phi radians around the given axis.
	 * 
	 * The returned rotation matrix can be used for three dimensional rotations in the euclidean space. If the given
	 * axis is not yet normalized, this method will normalize it.
	 * 
	 * @throws IllegalArgumentException if the given axis is the zero vector or has not exactly three components.
	 */
	public static MathMatrix get3DRotationMatrix(double phi, MathVector axis) {
		if (axis.norm() == 0.)
			throw new IllegalArgumentException("Axis for rotation matrix must be a non-zero vector.");
		if (axis.dim != 3)
			throw new IllegalArgumentException("Axis for rotation matrix must have exactly three components.");
		axis.normalize();
		
		// construct rotation matrix
		MathMatrix rotation = new MathMatrix(3, 3);
		double cosPhi = Math.cos(phi);
		double sinPhi = Math.sin(phi);
		double x = axis.get(0);
		double y = axis.get(1);
		double z = axis.get(2);
		rotation.set(0, 0, cosPhi + x*x*(1-cosPhi)  );
		rotation.set(0, 1, x*y*(1-cosPhi) - z*sinPhi);
		rotation.set(0, 2, x*z*(1-cosPhi) + y*sinPhi);
		rotation.set(1, 0, y*x*(1-cosPhi) + z*sinPhi);
		rotation.set(1, 1, cosPhi + y*y*(1-cosPhi)  );
		rotation.set(1, 2, y*z*(1-cosPhi) - x*sinPhi);
		rotation.set(2, 0, z*x*(1-cosPhi) - y*sinPhi);
		rotation.set(2, 1, z*y*(1-cosPhi) + x*sinPhi);
		rotation.set(2, 2, cosPhi + z*z*(1-cosPhi)  );
		
		return rotation;
	}
	
	/**
	 * @brief Returns a rotation matrix for rotations of phi radians around the given axis.
	 * 
	 * The returned rotation matrix can be used for three dimensional rotations using homogeneous coordinates. If the
	 * given axis is not yet normalized, this method will normalize it.
	 * 
	 * @throws IllegalArgumentException if the given axis is the zero vector or has not exactly three components.
	 */
	public static MathMatrix get4DRotationMatrix(double phi, MathVector axis) {
		if (axis.norm() == 0.)
			throw new IllegalArgumentException("Axis for rotation matrix must be a non-zero vector.");
		if (axis.dim != 3)
			throw new IllegalArgumentException("Axis for rotation matrix must have exactly three components.");
		axis.normalize();
		
		// construct rotation matrix
		MathMatrix rotation = new MathMatrix(4, 4);
		double cosPhi = Math.cos(phi);
		double sinPhi = Math.sin(phi);
		double x = axis.get(0);
		double y = axis.get(1);
		double z = axis.get(2);
		rotation.set(0, 0, cosPhi + x*x*(1-cosPhi)  );
		rotation.set(0, 1, x*y*(1-cosPhi) - z*sinPhi);
		rotation.set(0, 2, x*z*(1-cosPhi) + y*sinPhi);
		rotation.set(0, 3, 0.                       );
		rotation.set(1, 0, y*x*(1-cosPhi) + z*sinPhi);
		rotation.set(1, 1, cosPhi + y*y*(1-cosPhi)  );
		rotation.set(1, 2, y*z*(1-cosPhi) - x*sinPhi);
		rotation.set(1, 3, 0.                       );
		rotation.set(2, 0, z*x*(1-cosPhi) - y*sinPhi);
		rotation.set(2, 1, z*y*(1-cosPhi) + x*sinPhi);
		rotation.set(2, 2, cosPhi + z*z*(1-cosPhi)  );
		rotation.set(2, 3, 0.                       );
		rotation.set(3, 0, 0.                       );
		rotation.set(3, 1, 0.                       );
		rotation.set(3, 2, 0.                       );
		rotation.set(3, 3, 1.                       );
		
		return rotation;
	}
	
	/**
	 * @brief Returns a unit (n x n) matrix.
	 */
	public static MathMatrix getUnitMatrix(int n) {
		MathMatrix matrix = new MathMatrix(n, n);
		for (int i=0; i<n; i++)
			matrix.set(i, i, 1.);
		return matrix;
	}
	
	/**
	 * @brief Returns the product of the given matrices.
	 * @throws IllegalArgumentException if m1.m != m2.n (matrices can not be multiplied).
	 */
	public static MathMatrix mul(MathMatrix m1, MathMatrix m2) {
		if (m1.m != m2.n)
			throw new IllegalArgumentException("Given matrices can not be multiplied (dimensions missmatch).");
		MathMatrix result = new MathMatrix(m1.n, m2.m);
		for (int i=0; i<result.n; i++) {
			for (int j=0; j<result.m; j++) {
				double sum = 0.;
				for (int k=0; k<m1.m; k++)
					sum += m1.get(i,k)*m2.get(k,j);
				result.set(i, j, sum);
			}
		}
		return result;
	}
	
	/**
	 * @brief Returns the requested matrix - vector product.
	 * @throws IllegalArgumentException if
	 */
	public static MathVector mul(MathMatrix matrix, MathVector vector) {
		if (matrix.m != vector.dim)
			throw new IllegalArgumentException("Cannot compute matrix-vector product (dimensions missmatch).");
		MathVector result = new MathVector(matrix.n);
		for (int i=0; i<result.dim; i++) {
			double value = 0.;
			for (int j=0; j<vector.dim; j++)
				value += vector.get(j)*matrix.get(i, j);
			result.set(i, value);
		}
		return result;
	}
}
